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I. Introduction 


This study considers the developnent of a digital 4D automatic control 
law to capture and follow a steep glideslope (6°) under low visibility 
conditions and in turbulence, using the Microwave Landing System (MLS) 
under developnent by the FAA. The study of curved 4D flight paths leading 
to a steep final approach under low visibility conditions is part of the 
Terminal Configured Vehicle program (TCV), sponsored jointly by NASA and 
FAA. The goals of the TCV program include the reduction of aircraft noise 
in airport communities, the reduction of fuel consun 5 >tion, the reduction 
of the effects of adverse weather conditions on aircraft operations in 
air terminals, and the efficient use of airspace in congested terminal 
areas through the use of the Microwave Landing System. 

The major effect of the use of steep glideslopes is in the area of 
noise reduction. In comparison to the currently used 2.5° to 3° ILS 
^ideslopes, the 6° glideslope reduces the noise perceived on the ground 
due to its altitude profile and thrust level. At equal distances from 
the runway, the altitude of an aircraft following a 6° glideslope is 
almost twice the altitude of an aircraft following a 3° glideslope. 

Thus, the noise level heard on the ground is reduced due to the difference 
in altitude even when the same amount of noise is generated by both air- 
craft. A further reduction in noise is due to the fact that the aircraft 
flying the 6° glideslope generates less engine noise, so this steep 
glideslope requires a lower thrust setting than vjould be required by the 
same aircraft flying a 3° glideslope. This reduction in thrust level is 
of the order of 2:1 for the RSFS aircraft Of the TCV program. 



The reduction in thrust level associated with steep glideslopes 
also reduces the fuel consumed during the final approach. The ability 
to fly varying glideslope angles may also provide a method to avoid the 
vortex generated by large aircraft, by allowing smaller aircraft to fly 
different glideslopes to reduce the likelihood of such encounters ; 
however, further research in this area is necessary. In general, the 
ability to fly steep glideslopes provides a versatility that can be 
useful in efficient use of airspace in the terminal area. 

The guidance information necessary to fly steep glideslopes in low 
visibility instrument approaches can be obtained from the Microwave 
Landing System (MLS). The MLS is a ground-based guidance system which 
provides position information to aircraft inside its volumetric coverage. 
It consists of a DME providing range information, an azimuth antenna 
colocated with the DME providing the aircraft's azimuth angle relative 
to the runway up to ±60°, and an elevation antenna located at the glide- 
path intercept point but offset to the side of the runway providing the 
aircraft's elevation -angle up to 20°. A second elevation antenna located 
further down the runway to provide flare guidance is also under con- 
sideration. The MLS thus has a volumetric coverage, and provides guidance 
information that can be used for steep approaches and curved flight 
paths. The major characteristics of the MLS include high accuracy of 
position information, low sensitivity to adverse weather conditions and 
volume-trie coverage. 

With the high accuracy in position information provided by the MLS, 
it is of interest to investigate -the use of low accuracy accelerometer 


2 



data in place of niore sophisticated and costly systems such as inertial 
platforms in automatic landings under turbulent weather conditions. 

Thus, in this study, body-mounted accelerometers were used to provide 
acceleration information. This information was mixed with MLS data 
arriving at discrete instants of time and with air data in a constant 
gain Kalman filter to obtain optimal estimates of the aircraft velocity 
and sink rate as well as the wind velocities by filtering out the noise 
associated with the various sensors. The development of this filter 
for the longitudinal axis is given in Section III.B. The results 
obtained from a simulation of the filter are shown in Section V. 

In Section II, the aircraft’s equations of motion used in the 
simulation are described. A mathematical model describing the deviations 
of the aircraft’s longitudinal variables from their steady values on a 
6° glideslope is obtained. The effect of lags in thrust build-up and 
the effects of winds on the aircraft motion are included in this model. 
Using the Dryden spectrum, a dynamical model for the simulation of 
wind gusts is developed, then steady winds are added to this model. The 
models are expressed in state variable form which is more suitable to 
the use of modern estimation and control techniques. 

In Section III, a mthematical model for the noises in the various 
sensors is developed. Then a non-linear pre-processor is used to 
transfoirm these measurements into a form more suitable for filtering 
purposes. In Section III.B, the development of the filter is described 
and some aspects of its implementation are discussed. 
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In Section IV, a digital automtic control law to capture and follow 
a 6° glideslope is developed for the longitudinal axis. The control 
law uses the aircraft variables as well as the wind estimates to decrease 
the aircraft’s deviations from the glidepath. 

Section V describes the results obtained from a simulation of the 
aircraft, winds, sensor errors, the filter and the control law. 

It is a pleasure to acknowledge Dr. Thomas M. Walsh for his en- 
couragement of the concepts presented in this study. 
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II. Modelling of Aircar^t Dyrianiics and Winds 

The general equations of motion for aircraft are complex nonlinear 
differential equations and can be found in various texts on aircraft 
dynamics [4], [53, [63. As this study considers the glideslope capture 
and glideslope. phases of the final approach, however, several sinplifying 
assumptions can be made to reduce the complexity of these equations and 
make them more amenable to analytic manipulations without appreciable 
degradation in their validity [4, p. 2303, [5, pp. 254-2653. As the 
equations of motion are used extensively in the study and the computer 
simulation, the specific equations used will be described here. 

A. Aircraft dynamics with wind disturbances 

The phase of flight considered in this study is glideslope capture 
followed by a steep glideslope up to 6 degrees flown at a constant 
airspeed of 120 knots except for small fluctuations. In these phases of 
flight (Fig. 1) the aircraft is aligned with the runway, has a zero 
or very smll yaw angle with respect to the runway as well as a zero 
bank angle except for the case of a significant cross-wind requiring a 
"crab" maneuver. The control activity for the lateral motion is aimed at 
keeping the aircraft aligned with the runway, with level wings. Hence, 
all the lateral variables, i.e. yaw, roll, their rates and the sideslip 
angle, have very snail values except for crab maneuvers. Similarly, 
among the longitudinal variables, the pitch angle is small during these 
phases of fligjnt, usually within 6° to -4° for a 6° capture. 

Under these conditions, the nonlinear equations of motion can be 
linearized about the steady state flight condition of the glideslope 
using well-]<ndwn methods [43, [53, [63. The deviations from the steady 
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flight condition can be described by linear differential equations 
which are simpler to use in analytical manipulations and are more 
suitable for the application of modem control tdieory principles. Under 
the conditions stated above (these conditions will be restated more 
precisely below), the equations of motion of the aircraft can be 
ej<pressed as [2, p. 2.32]: 


m(u - V r - R^v + W^q + Q^u)) = -mg cos0^9 + f^ ’ 

X X 


( 1 ) 


m(v +Ur + Ru-Wp-Po)) = -mg sin$ sin0 0 
o o o o ^ o o 


+ mg cos$^cos0^(() + f^ ’ 

y y 


( 2 ) 


m(u) - Uq-Qu + Vp + Pv) = -mg cos$ sin0 0 
o o o ^ o o 


- mg sin$ COS0 (f> + f . + f™ , (3) 

o o A T 


I p-I r-I (Pq + Qp) + (I -I )(Rq + Qr)=£.. + 
xx^ xz xz o^ ^o^ zz yy o o A T’ 

(4) 


I q + (I - I )(P r + R p) + 21 (P p - R r) = m. + nv-, 
yy^ XX zz o o^ xz o^ o AT’ 


(5) 


I f-I p+(I -I )(Pq + Qp)+I (Qr + Rq)=n. + n^, 
zz xz^ yy XX o^ ^o^ xz o o A T 


( 6 ) 
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where 


= steady state inertial speed in the x direction 
= steady state inertial speed in the y direction 
s steady state inertial speed in the z direction 

u = perturbation in the inertial speed in the x direction 

V = perturbation in the inertial speed in the y direction 

CD = perturbation in the inertial speed in the z direction 

= steady state roll rate 
steady state pitch rate 
= steady state yaw rate 
p = perturbation in roll rate 
q = perturbation in pitch rate 
r = perturbation in yaw rate 
i] steady state roll angle 
0^ = steady state pitch angle 
B steady state yaw angle 
<|) = perturbation, in roll angle 
9 = perturbation in pitch angle 
ij) = perturbation in yaw angle 

= perturbation in net aerodynamic force along the x direction 

^Ay - perturbation in net aerodynamic force along the y direction 

f/^ = perturbation in net aerodynamic force along the z direction 

f^^ = perturbation in thrust along the x direction 
^Ty “ perturbation in thrust along the y direction 
frp^ = perturbation in thrust along the x direction 
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pertxirbation in rolling inoment due to aerodynamic foroes 
perturbation in pitching moment due to aerodynamic forces 
perturbation in yawing moment due to aerodynamic forces 
perturbation in rolling moment due to thrust 
perturbation in pitching moment due to thrust 
perturbation in yawing moment due to thrust 

These equations are valid for any set of ri^t-handed rectangular 
body-fixed axes, i.e., right-handed reference frame, fixed to the body of 
the aircraft with the origin located at the aircraft's center of mass. 
Figure 2 shows the sign conventions and the vectors pictorially for the 
vertical plane. The assunptions and approximations used in arriving 
at equations (1) - (6) are given below. 

1. The earth is assumed to be flat and fixed in an inertial reference 
frame; 

2. The perturbations in the angles are siihII, so that 

cos 0 = cos<}> = 1 (7) 

sin Q - Q, sin<J) - (8) 

3. The second order terms in the perturbations are negligible 
relative to the first order terms in the perturbations; 

4. The aircraft is a rigid body. 

In this form, equations (1) - (6) are coupled; however, if the 
steady state fli^t condition is taken to be the glideslope, these 


= 

rVp = 
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equations can be decoupled and sirrplified. This steady state fli^t is 
described by: 


$=P=Q=R=V=0 (9) 

0 0^000 

Hence, for this steady fli^t condition, the equations of notion sinplify 
to: 


m(u + W q) 
o^ 


-mg cos$^0 + f^ + f^ 


m(u - U q) 
o^ 


-mg sin$^e + f^ + f^ 
z z 


V = -A * 


m(v + Ur-Wp) = mg cos$ * + f . + f_ 

o o^ ^ o^ A T 

y y 


I p-I r = A. +t 
xx^ xz AT 


( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 


I r-I p = n. + ru, (15) 

zz xz^ A X 

Equations (10) - (12) contain only longitudinal variables, whereas 
Eqns. (13)- (15) contain only lateral variables so that the equations are 
now decoupled. As this study is concerned with glideslope capture and 
glideslope tracking, only the longitudinal equations of motion (i.e. , 
Eqns. (10) - (12)) will be considered in the following. 
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Throu^out this study, mainly tiiree sets of coordinate axes will be 
used: the earth fixed axis, the body axis and the stability axes. The 

earth-fixed coordinate frame (x^, y^, z^) has its origin fixed on a 
specified point on the runway at which the aircraft is going to land; 
the x^ axis is along the runway, the direction in which the aircraft 
will land being positive; i.e., at touchdown the A/C will have a positive 
veolcity along x^. The z^ axis is vertical positive downwards; and y^ 
is perpendicular to both x^ and z^ with the positive end in the direction 
to make the coordinate frame right-handed. As the earth is assumed to 
be stationary with respect to an inertial frame, (x^, y^, z^) is itself 
an inertial frame. 

Two body-fixed axes with their origin fixed at the center of nass 

of the A/C are also used: the body (Xj^, y^, z^) and stability (x^, y^, z^) 

axes. The x^ axis is in the A/C’s plane of symmetry and is taken 

to be along the fuselage reference line of the A/C, positive towards the 

nose; the y^ axis is positive towards the ri^t wing, and is positive 

downwards; this reference frame will be referred to as the body axes. 

The stability axes (x , y , z ) can be obtained from the body axes by a 

s s s 

rotation of about the y^^ axis such that when the A/C is in its steady 
state fli^t condition, its velocity vector is along the positive x 

s 

axis. The equations of motion for the aircraft (10 - 12) will be written 
in the stability axes. 

In the above equations, the term f^ represents the total algebraic 

X 

change in the value of the aerodynamic force along the x axis due to 
changes in the values of the aerodynamic and control surface variables; 
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the changes in the forces and variables are referenced to the steady 
state values of the corresponding forces and variables. The remining 
terms in equations (10 - 12) are defined similarly. For the longitudinal 
eqns. of motion, the aerodynamic forces and moments involved are the 
lift and drag forces and the pitching moment; the effects of thrust are 
described separately by the terms with subscript T. These forces and 
moments can be expressed as follows: 


L = Cj^(u , a, a, q, 6e, 6s)qS (16) 

D = Cj^(u, ot, a, q, 6e, 6s)qS (17) 

T = C^(u, a, a, q)qS (18) 

M = Cj^(u, a, a, q, 6e, 6s, 6T)qSc (19) 



where L, D, T and M represent the lift, drag, thrust forces and the 
pitching moment respectively, q is the dynamic pressure, V the airspeed, 
u the perturbation of V *s conponent along the x axis, a the perturbation 
in the angle of attack, a the time rate of change of a, q the pitch rate 
of the A/C relative to the atmosphere, 6e the perturbation in 1die elevator 
surface deflection, 6^ the perturbation in the stabilizer surface de- 
flection, 6T the perturbation in the thrust force, S the wing area and 
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c the mean aerodynamic chord. Note that the aerodynamic forces and 
moments depend on the motion of the aircraft relative to tiie atmosphere 
such as u, a, etc. , rather than the inertial motion variables u, to, etc. 
Thus the effect of winds is automatically included into the equations of 
motion. If the pertijrbation in the angle of attack of the inertial 
velocity, say a, is defined as 


A , — 1 

a = tan 

o 


w 

+ u 



u « U , u « U , 
o o’ 


( 21 ) 


then 


u 

u = u + u,u' = j 7 -=u' +u* (22) 

w’ - Uq w 

a = a + a , (23) 

w 

ot = a + a , (24) 

w 

q = q + (25) 


where the subscript w denotes the component due to winds. 


Using Figure 2, it is seen that the total aerodynamic forces 


along the x and z axes are 
s s 


= Lsinot - Dcosa (26) 

s 


F. = -Lcosa + Dsina (27) 

Az„ . - 
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The foipces due to thrust are 


*■ 1 X 3 = 


( 28 ) 


^Tz = -Tsin<»o *T>' 

s 


(29) 


Fran these equations, we can obtain expressions for the perturbed 
force and noments by using a Taylor series where second and hi^er 
order terms are neglected, hence. 


3F 


Ax 


3F 


Ax 


3F 


Ax 


3F 


Ax 


3F 


Ax 


3F 


Ax 


Ax 


3u - 


u + 


3a 


a + 


3(5 


a + 




36e 


6e + 


36s 

(30) 


6s 


where the partial derivatives are evaluated at the steady state condition. 

Similarly, f^ , f^ , , m^ and nv^ can be expressed by a truncated 

s s s 

Taylor series of the same form. The partial derivatives in these 
expressions can be expressed in terms of the partial derivatives of 
the lift, drag, thrust and pitching moment coefficients, i.e. , the 
stability derivatives. The equations of motion thus obtained are given 
below. 


1“ = -mg =ose^6 + %S((-Cj^, + + C^.H- 2C^)u' 


(31) 


^ ^^Lo ^6e^® " ^D6s^®.'*‘ ‘^Tx6T^'^^ 
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mCw - U^q) = -rag sinG^e + ^S{-(C^,+ 2C^)u» - (C^ + C^)a 

(32) 

- C, .a - C, qH-C. 6e - C, 6s + CL,^ 6T} 

V = ^ * 'Sb " * <k“ 

(33) 

" V " =MSe*® * >^61^1}, 

where has been taJcen to be zero. After some manipiilation, these 
equations can be expressed in state variable form. However, we shall 
first derive some further equations and include those as well in the 
state variable raodel. 

First note that the flight path angle, y, shown in Figure 2 <3an 
be expressed as 

Y=0 +0-(a +a)=(0 -a)+(0-a)=Y +(9-a) 

' o o o o o 

(34) 

Now, let (x,z) be the coordinates of the aircraft's center of mass 
along the earth-fixed coordinate frame x^ and z^. Then the ground speed, 

X, and the sink rate, z, are given by 

X = V.COSY = V.cos(y +0-a) (35) 

X X o 

= V.cosacos(Y + 0) + V.sinasin(Y + 0) = Ucos(y + 0) + Wsxn(Y + 0) 

X o X 'o o o 

(36) 

z Q -V.sinY = -Usxn(Y + 0) + VhosCy^ + 0) 

^ ° 14 


( 37 ) 


since the stability axis is chosen to cororespond to the glideslope, 
the ground speed, x , and the sink rate, z^, of the aircraft in the 
steady state condition are given by 


X = U COSY > 
o o 'o 


(38) 


z = -U sinv 5 (39) 

o o 'o 


as in this steady state condition the speed along the z^ axis (i.e., W^) 
is zero. Hence, the perturbation in the ground speed and sink rate are 


x-x =(x-x) = Ucos(y + 0) + Wsin(y + 0) - U cosy (40) 
o o o o o o 


z-z n(z-z) = -Usin(Y + 0) + Wcos(y ‘ + 0) + U siny (41) 
o o o ‘o o o 


Substituting the equations 


U = U +u, W = W + u) = ( 0 , 

o ’ -o ’ 


(42) 


into (40) and (41), 


(x - X ) 
o 

U 


[cos(y + 0) - cosy 3 + u’cos(y + 0) + — sin(y + 0) 
'o 'o 'o U o 

o 

(43) 


z - siny^0 + cosy^u’ + siny^a; 


(44) 
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CsuiCy + 0 ) - sinv ] .- u'suiCv + 6) + n~ cos(v + 9) 

‘o o o U o 

o 

(45) 

- -cosy^S - siny^u’ + cosy^a. (46) 

Thus the perturbations in ground speed and sink rate can be 
expressed as linear functions of 0, u’ and a by (44) and (46) as long as 
the perturbations in the pitch angle, 0 , are sinall. Note that the 
perturbations (x(t) - x^(t)) and (z(t) - z^(t)) in ground distance and 
altitude correspond to the variable time, i.e.x(t) - x^(t), is the 
perturbation or error in the x position of the aircraft relative to where 
it should have been at that time. Hence, if these errors are zero, 
the aircraft not only keeps an average inertial speed of U^, but also 
has to be at specified positions at the appropriate time; i.e. x - x^, 
z - z^ represent 4D (four dimensional) errors which include time as a 
variable. 

To include the effects of the servo responses of the actuators, we 
shall also model the thrust and stabilizer motions dynamically by linear 
equations . Thus , 

6 T = -.56T + .2986th (47) 

6 th = U 3 , 6 *s = U 2 (48) 

where 6 T, 6 th and 6 s are perturbations of thrust, throttle and 

stabilizer from their steady state value, respectively. The elevator 
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is not modeled in this fashion as its motion is relatively fast, i.e. 
its time constant is much lower than the others. Thus the thrust model 
takes the "spool up" time of the engine into account, at least linearly. 

Now, the equations of motion of the aircraft, the 4D perturbations in 
the ground speed and sink rate and the effects of the actuator responses 
can be combined into a state variable model. Let the state vector x be 
defined as 


X - X z - z 

x' = (0 u’ a q — jj— ^ -- 6T 6th 6s) (49) 

o o 

Then equations (31) - (33), (44), (46) - (48) can be combined into the 
vector equation 

X = Ax + Bu +Dw (50) 

vhere u' = (6e 6s 6th), w = (u^ q^) and A , B and D correspond 

to the appropriate coefficients in the original equations and are 
given below 


0 

0 

0 

1 

0 

0 

0 

0 

0 

a.21 

^22 

^23 

0 

0 

0 

^27 

0 

^29 

^31 

^32 

^33 

^34 

0 

0 

^37 

0 

^3 9 

am 

a42 

^43 

a44 

0 

0 

a47 

0 

^49 

-siny 

o 

cosy 

o 

siny 

O 

0 

0 

0 

0 

0 

0 

-cosy^ 

-siny 

o 

cosy 

o 

0 

0 

O' 

0 

0 

0 

0 

0 

0 

0 

0 

0 

ay 7 

^78 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

_ 0 

0 

0 

0 

0 

0 

0 

0 

0 
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0 

0 

0 " 


0 

0 

0 *" 

b2i 

0 

0 


^22 

3-23 

0 

t>31 

0 

0 


332 

^33 

d33, 

Wl 

0 

0 


a42 

^43 

<^43 

0 

0 

0 

, D = 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

1 


0 

0 

0 

0 

1 

0 - 


0 

0 

0 _ 


The expressions for the elements of the rfetrices A, B and D are 
given in the Appendix. Hence, the aircraft equations of motion along 
with 4D error equations in ground speed and sink rate have been modeled 
by a linear state variable model. It sould be noted, however, that the 
equations have been obtained by linearization of 1±ie nonlinear equations 
about the steady fli^t condition of a glideslope with angle and a 
constant airspeed U^; hence, these equations provide a realistic model 
of the aircraft motion provided that the deviations from this steady 
fli^t condition are small. 

B. Wind ?fodeling 

To complete the aircraft model given by (50), the wind vector w has 

to be specified. The components of this vector consist of u* , the 

w 

nortiHlized wind velocity in the -x^ direction; a^. the part of the angle 

of attack due to winds, and q^, the rotation of the atmosphere about the 

y^ axis. The and conponents are modeled as consisting of a gust 

conponent with zero average value and a steady component (i.e. the 

average value u' and a ). Hence, 
w w 
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u* = u’ + u* = (u + u )/V 
w g s g s a 


( 51 ) 


w 


a + a = -ill! + 0) )/V 
g s g s a 


(52) 


where is the airspeed. 

The gust con^xDnents can be modeled using the well-known Dryden 
spectrum [4]. This method consists of using spectral factorization methods 
to obtain a dynamical system which generates a random process having the 
specified power spectral density when driven by white noise. Let u (t) 
be the gust velocity with respect to earth at time t; then the covariance 
function of the random process u (t) is defined as 


R (t) = E(u (t)u (t + t)), 
u_ g g 


(53) 


where ’^g('t) is assumed to be a wide-sense stationary random process 
with zero mean and E denotes the statistical expectation operator 
[7], [8]. The power spectral density of this process is then defined 

i': 

as the Fourier transform of its covariance function R (t): 

u_ 


S (o)) = / R (t)e‘^“'^dt; j = /^T; 
g g 


(54) 


then, the variance of u (t) or the power of the random process is given 
by 


* Definitions of the Fourier transform differing from (54) by the 
factors l/2ir, 1/v^^, 2/ir are sometimes used in the literature, we 
shall use the definition given above. 
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R (0) = ^ f S ((o)do). 
u 2ir u 

g g 


(55) 


The Dryden spectra describe the statistical behavior of the wind 
gust velocities in the aircraft body coordinates by specifying their 
power spectral densities in terms of the spatial frequency [9]. 


2 

2a L 

q (-o') - u u 
u " 1 + (L 


a L 1 + 3(L 
o ('O'* - “ “ “ 

UTTLliKF 

a 03 

Ji2v2 

•s (J 2 ) = s (J 2 ) , 

1 a)^ “g 


(56) 


(57) 


(58) 


vtiere b is the wing span, and are the scales of turbulence, and 
is the spatial frequency related to the temporal frequency a* by 


= (o/V 


(59) 


The wind gust velocities along the aircraft's body axes will be 

denoted by the subscript b. It is assumed that u is uncorrelated with 

§b 

both u) and q ; but w and q are correlated, since q is due to 

Sb §b Sb % Sb 

the variation of a along the aircraft’s body. Using Taylor’s hypothesis 
Sb 

of a "frozen field," [10], these spectra can be expressed in terms of tlie 
temporal frequency w by 


S^(o)) = ^ i = u, ( 0 , 

0. 0. 


(60) 
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A mathematical model of the gusts can be obtained using spectral 
factorization [11], [12], [13]. 

1. Model of u^ - Using (56), (60) and (51), the spectrum of u’ 

^ g] 

can be found to be 


S^,((o) 


2L , 

u u 1 

L 2 2 

^ 1 + (^) 0) 
''a 


(61) 


2L 
u u 

a 


(62) 


T • U T . • u 

1 - 1 + 

a a 


This corresponds to a system with transfer fijnction G^(s) driven by 


a white noise with power 2L / V^(ft^/sec^)/Hz. 

^ u u a 


G (s) 
u 


L 

1 + ^s 

a 


(63) 


2. Model of agj^ - Using (57) and (60) the spectrum of can be 


foimd to be 


1 + 3 (^) 0 )^ 

= °ivl r 17^2 


"" L 2 

1 + Cy— ; oi 
_ ^ 


(64) 


0) 

~W 


1 - 


3 . 


1 + j/3 ^ 0) 
''a 


3. • ' 3 ' 


( 65 ) 
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From equation (65), it is seen Itiat a system with transfer function 
G^(s), given below, driven by white noise with power c^L^/V^(ft2/sec2)/Hz 
generates the spectrum in (64). 

lt/3^s 

a 


G (s ) = 
a 




1 + 2 


( 66 ) 


3. Model of qgj^ - Using (58) and (60), the power spectral density 


of q can be found to be 

% 


S ((d) 

q 


♦ S (to) 




(67) 


1 + ( 


irV 


G (-jo)) • — G (juj) . 

T .4b n,-4b a-* 

a a 


( 68 ) 


The cross-correlation of q and a is specified by their cross- 


Si 


Sb 


spectral density 


S (to) = 
qa 


30) 


~)?K S (to) . 

• 4b . . a 


(69) 


1 


From ( 68 ) and (69), it can be seen that if a (t) is input into 

% 

the filter 


~ \b 
^ 1 + - 77 - s 

nV 


( 70 ) 
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then, both the power spectral density and the cross-spectral density 
requirements will be met. A block diagram for generating the wind gust 
velocities is given in Figure 3. The inputs Wj(t) and W 2 (t) are 
uncorrelated gaijssian white noise processes. It should be noted that to 
obtain a stationary process, an initial period for the settling of the 
transients due to unmatched initial conditions need be allowed or the 
initial conditions so chosen as to obtain a stationary process 
immediately. Otherwise, the process would be non-stationary for a time 
corresponding to a few time constants. 

4. ffodels of u’ and a - The steady winds can also be modeled by 
s s 

differential equations, i.e. , by setting the derivative of the variable 
equal to zero, and the initial condition equal to the value of the steady 
wind. To allow for slow variations a forcing function with a small 
magnitude can be added. Thus, the steady components with respect to 
the earth-fixed axes are modeled as 

u ( 0 ) CO ( 0 ) 

u' □ W 3 , = Wc, ; u!(0) = -g , a^(0) = -g (71) 

o o 

vhere u ( 0 ) and to ( 0 ) are the steady wind values, and W 3 , w^ are in- 
s s 

dependent white noise processes with very small power. To, obtain the 
total winds u^ and a^, the gust and steady <x>mponents are first 
transformed into the stability axes and then added as shown in (51), (52). 

Thus the wind vector w in (50) can be expressed as the output of a 
linear system. The transfer functions G (s), G (s) and G (s) can be 
expressed in differential equation form, and then these equations can 
be put into state variable form. 
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(72) 


W = V + w = C^W, 


where A , B and C are given by 
w w w 


A = 
w 


0 

v2 

a 

~I7 

0) 

ttV 

3 . 

4b 

0 

0 

0 


2V_ 

LoT 

0 

0 

0 

0 


irV 

a 

4b 

0 

0 

0 


V. 


u 


0 

0 


0 

0 

0 

0 


0 

0 

0 

0 


B 


w 


V 

“ V 2 

(1 - 42 )(|^) 

OJ 

0 


0 

0 

0 


0 

0 

0 


V 


01 

0 


0 

0 

0 

0 

1 

0 


0 

0 

0 

0 

0 

1 


■c = 

w 


-sina 0 

o 


-cosa 


ttV 

a 

4b 


0 -cosa -cos(o +0) sin(a +6) 
o o o 


0 

ttV 

a 

4b 


sina -sin(a +0) -cos(a_ + 0) 
o o o 
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Thus, equation (50) describes the aiixrraft's notion due to changes 
in the control variables and wind conditions, while (72) describes the 
statistical properties of wind gust velocities and steady winds. Hence, 
the aircraft's equations of notion have been e^ressed by a linear d 3 ?nand.cal 
system of equations in state variable form. 

C. Discretization of the Equations of Motion 

In this section, the equation of notion given in (50) and the wind 
equations given in (72) are discretized; i.e. , the differential equations 
(50), (72) are replaced by difference equations which update the state 
variables from one sampling instant to the next. There are two major 
reasons that lead to the discretization of the equations of notion. 

The first is due to the simulation of the aircraft's notion on a digital 
con 5 )uter. Due to the nature of digital computers, the integration of the 
equations of notion has to be performed in a discrete manner. The 
second, and more important, r^son is inherent in the operation of 
some of the components of the system. The aircraft's position is obtained 
from the Microwave Landing System (MLS) which provides this data at discrete 
intervals of time, rather than continuously. Furthermore, the aircraft's 
control system includes a digital computer to perform the operations 
required for the implementation of the control law. Thus, thfe control 
commands at the computer output are, of necessity, discrete. These 
considerations lead to the discretization of the equations of motion. 

Let tj^ = !kT be the times at which the aircraft's state has to 
be known; then it is required that the state be updated from tj^ to tj^+2* 

This can be done by integrating equation (50). The result is given by 
[14]: 
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\+l 

^ ‘*’^''^k+l " s)Bu(s)ds 
''^k 

Vl 

+ / - s)EW(s)ds, 


where (j>(t) is the transition matrix of A; i.e. <f)(t) is the natrix 
A*t 

exponential e . Now, since the control law is digital, we shall assume 
that the control command does not vary over one sampling period, i.e.. 


^ ^ ^ \+i- 


Using (74) and the change of variable 


(73) can be expressed as 


[ T -IT 

J (j)(T - T)dTBj + / <j>(T - t)Dw(tj^ + T)dx. 


Now, the wind equations (72) can be integrated in a similar manner 


to obtain: 


"k+i = «k’ \ = ‘^w“k’ 
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where Wy. and Wj^ represent w(tj^) and W(t^), 3?espectively, and^ 


(j) = e ^ , C = / (J) (T - t)B 5(t + x)dT. 

k i ^w w T< 


It can be shown that { 5 ^^} is a white noise sequence [15], with covariance 


R = E(C, C> = / <f> (T - t)B FB»(},'(T - x)dT, 

C k k ■' ^w w w^w ’ 

o 


where F6(t - s) is the covariance of the continuoiis white noise process 
?(t). Now, note that 




w'"w'^'"k i ^w' 
o 


After substituting (79) into (75) and some manipiilation 


Vi ^ *\*^'\* 


where 5 ^ represents x(tj,), and 


( 1 ) = e , r =[/ 4>(T - s)ds]B,r = / <}.(T - s)DC (}> (s)ds, (81) 

' W ' w w 

o o 


T T 


o s 


\ = / / (J>(T - x)DC^(j)^(x - s)dx £;(tj^ + s 


* The matrices A , B and C have been evaluated at the steady fli^t 
condition for these derivations. 
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It should be noted that {t), } is also a white noise sequence, and 

ic 

v^ch includes those effects of the winds which are not correlated to 
Thus, equations (80) and (76) are the discrete versions of the 
aircraft equations of motion (50) and wind equations (72), respectively. 
These equations can be programmed on a digital computer to simulate the 
motion of the aircraft under various wind conditions, and in the develop- 
ment of a digital control law for glideslope capture. 



III. Development of the Filter Equations 

The aircraft position is obtained from the Microwave Landing 

System (MLS). This data is obtained at discrete intervals of time and 

witii hi^ accuracy. The discrete chamcter of the data nakes it 

suitable for digital processing. Since the data has a higb accuracy, 

it is of interest to investigate the possibility of deriving the 

signal parameters used by the guidance and control system without the 

use of costly inertial navigation systems. This can be done by 

processing MLS data and other sensor outputs throxagh a filter which 

estimates the desired parameters. Hence, in this section, the equations 

defining the filter processing will be developed using MLS data, air 

data and body mounted accelerometers without using inertial platform 

data. A discrete-time Kalman filter with constant coefficients to i?educe 

the amount of on-board computation will be used for filtering purposes. 

A. Development of Measurement Mcdels 

To describe the characteristics of the data obtained from various 

sensors, sinple models which describe sensor errors statistically will 

be developed. Models for these measurements including the associated 

errors are required both for simulation and filtering developments. As 

a Kalman filter will be used for processing, the errors in the measurements 

will be described as additive for filter design purposes; however, when 

this is unrealistic a different model will be used for the simulation 

of the errors which are actually input to the filter. It will be 

assumed that the bias eirror in the measurements is either negligible or 

has been subtracted from the data; otherwise, a pre-filter should be 

included for that purpose. As this study considers the longitudinal 
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eqmtions of motion, the measurements considered here are those v^ch 
affect Idle longitudinal variables. The measurements considered and the 
associated source of the measurements are listed below. 


Yj = pitch angle (gyro) 

Y 2 = pitch rate (gyro) 

Y 3 = slant range (MLS) 

Y^ = elevation angle (MLS) 

Yg = altitude (barometric) 

Yg □ sink rate (baronetric) 

Yy = acceleration along z (body-mounted accelerometer) 

o 

Yg □ true airspeed (air data computer) 

Yg = acceleration along (body-mounted accelerometer) 

The sensor models used for the simulation of the measurements 
given below. 

Yj = 0^ + 0 + VI = 0^ + xi + Vj 


Y 2 = 6 - V 2 = X 4 + Vj 


2 1/2 

Yo o C(U X5 + X - X ) + (U Xg + z - z )2] + V, 

^ o^oa o^oa 3 


'U Xg + z - z 
Y, □ tan-i ' ° ■ ■ 


Vs 


z - z \ 

o e I 

+ X I 

o / 


+ V 4 


are 

(83a) 

(83b) 

(83c) 

(83d) 
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Ys = -(U^xg + + V5) 


(83e) 


Yg = -CU^ie + (83f) 

1 + X 2 . 

Y7 = U [xotanXii + n — Xo 

' O ^ COS'^X3 ^ 

/ 1 + u’ \ 

Yg = U ( ^ ) Vg 

° o\ cosa / ° 

Yg = U^[x2 + (1 + X 2 )xijtanx 3 + Xgsinxg - x^sinx^] + Vg (83i) 

where x^ and are the 4D coordinates of the desired glideslope at a 

given time, x^ for i > 1 is the i^ component of the state vector x 

given in (50), is the vertical coordinate of the elevation antenna, 

X , z are the coordinates of the azimuth antenna and V. is the noise 
a. a. X 

introduced by the sensor. The expressions for the accelerations in 

Yy and Yg are obtained by writing the inertial velocities along the 

x^ and z^ axes in terms of the state variables, differentiating with 

respect to time and then transforming these accelerations to the stability 

axes, X , z . The earth-fixed coor>dinate system is referenced with 
s s 

respect to a point on runway centerline corresponding to the position of 
the elevation 1 antenna. The values for the standard deviations of the 
sensor noises were chosen to reflect current instrumentation standards 
[16], [17]; these are shown in Table 1. 

To use these measurements as input to the Kalman filter it is 
necessary to express them as linear ccaribinations of the aircraft state 


- X4(l + X2 + COSXj - COSX3)] + Vy 


(83g) 


(83h) 
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variables x and wind variables W with additive noise. To achieve 
this a pre-filter processor is used. This processor is nonlinear and 
consists of a general coordinate transfoimation to obtain the variables 
suitable for filtering. The equations describing the processing performed 
are given below. 


r = [Y^ - 2 L Y cos(Y + n)]^'^^5n = tan“^ ( 84 ) 

636034 ' 

Yl = Yi - 0^ = Xi + vi (85a) 

Y2 = Y2 = X4 + V2 ( 85 b) 

Y3 = 0^ [r^cosY4 + x^] = X5 + V3 ( 85 c) 

o 

74 = ^ [r^sinY4 “ + Zq] = xe + V4 ( 85 d) 

o 

75 = if [z^ - Y 5 ] = xg + V 5 (85e) 

o 

76 = jf - Yg] = e'Ax + V6. (85f) 

Uq o g 

77 = ^ - EgBu = e’(Ax + DC^W) + V7 ( 85 g) 

o 

78 = ^- 1 = X2 - eJC^W + Vg (86h) 

o 

79 = ^ = e^^Ax + DC^W) + vg ( 85 i) 

o 
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viiere h is the height of the elevation antenna above the azimuth antenna, 
ea 

X is the distance of the Az. antenna from the origin, L is given by 
a e 

(x /cosn)» f is the range fcom the elevation antenna, e. is a column 
a e 3 

vector equal to the colimin of the identity matrix and is a pseudo 
noise representing the additive noise in the processed measurement vector 
y. The first equality in (85) gives the processing done to obtain y 
fran the total measurements Y, while the second equality (i.e. the 
right hand side of (85)) provides the measurement model which is used 
in the development of the filter. This model can be written in matrix 
form as 


= yjc = ^ ^86) 

where and C 2 correspond to the appropriate coefficients given in 
the second equality of (85). Thus, a nathematical model to simulate 
the noisy measurements Y, a nonlinear processor and a linear measuranent 
model with additive noise for the development of a Kaljien filter are 
obtained. 

B. Developnent of Filter Equations 

The output of the pre-filter processor, y , given in (86) can be 

jC 

used as the input to the filter. In turn, the filter reduces the noise 
introduced by the sensors by optimally weighing new data versus previous 
estimates and also generates estimates of variables which are not 
directly measured using the equations by which the dynamical system is 
governed. The variables which are not dicectly measured include the 
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angle of attack and the wind velocities; thus the filter output provides 
optimal estimates of all the state variables and the wind velocities. 

The control system uses these estimates to calculate the control surface 
settings thus closing the loop around the aircraft viewed as a d 3 ?naniical 
system. 

Since estimates of the wind velocities as well as the aircraft state 
variables are to be obtained, it is desirable to combine the aircraft 
equations of motion (80), and the wind model (76) into a single system 
of equations. Thus, let 





In this notation (80), (76) and (86) can be written as 


(87a) 


(87b) 


* 'k* = “k * \ 


( 88 ) 


Now, it is desired to obtain estimates, Xj^, of using the 
available measurements y, . A Kalnan filter for the system given in (88) 
would provide optimal estimates of X^. Note that, whatever processing 
is done by the filter, it will require a certain amount of time on a 
digital computer; thus the estimates will be available for use by the 
control system only after the filtering computations have been 


performed. It is known that such delays in the control loop can lead 
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to instability [18] in digital control systems if they are not ccsnpensated. 
Thus it is desirable to obtain estimtes with nujumum delay. This can 
be accomplished by predicting from the past measurements 
^kl’ ^k2’ using the prosent measurement yj^. In this 

way,- at time "the measurements y^^ are obtained and the conputations 

necessary to obtain are initiated. If these computations can be 
performed in less than one sampling period, T, then Xj^ is available at 
time tj^. Hence, a one-step predictor algorithm is best suited for this 
application. The optinrol prediction algorithms for the system given 
by (88) are the well-known Kalman equations [19], [20], [21]. 


\+l ~ ^^k ^^^k 


(89a) 


-1 


^k+1 ” aPj^C*[CP}jC’ + R|^] CPj^a’ + Qj^ j 


(89b) 


Gk = aP3^C«[CP^C* I^] 


-1 


(89c) 


where P^^ is the covariance of the error, ^ "the measurement 

and state noise covariances are given by 


E(?k ej) = , E(Vj^ vp = 


(90) 


Note that (89a) is a predictor algorithm since it uses Uj^ which is 

y«i 

available at time tj^ to compute which is not used until (89b) 

propagates the error covariance and (89c) gives the filter gain matrix. 
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It is important to note that to update the estimates recursively as 
shown in (89a) only the gain matrix Gj^- is needed. 

It is well-Jcnown that if the noise covariances and do not 
depend on time, i.e. , 

o Q, \ ° Rj all k (91) 

then the solution of the matrix equation of Riccati type given by 
(89b) converges to a steady state error covariance, say P, for a broad 
class of dynamical systems (a,3,C). The observability of (a,C) is a 
sufficient but not necessary condition for the convergence of P, ; a 
necessary and sufficient condition is given in [21]. Thus, in the steady 
state situation, the gain also converges to a corresponding value, 
say G, given by (89c). Hence, the filter (89a) becomes a time- invariant 
system. The implications of this convergence to the implementation of 
the filter are that the error covariance equation (89b) can be solved 
off-line and the steady state gain G can be computed prior to fli^t; 
so that the only computations which need to be performed in real time, 
i.e. during fli^t, are the operations described by (89a). As the 
propagation of the error covariance generally requires a very large 
number of computations compared to the update of the estimates, the 
use of the steady state gain reduces the number of operations that need 
to be performed on-line considerably. The reduction in coiputation time 
thus gained nay allow the use of a sophisticated model for the aircraft 
dynamics, thus increasing the accuracy of the estinates. The number 
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of computations required can further be decreased by taJdng advantage 
of the special properties of the matrices involved. However, for the 
most general case where no special properties are present, the 
nxjmber of multiply-add operations required to update the estimates fran 
one sample to the next is roughly given by n(n + p + 2m) , where n is 
the number of variables vhich are estimated (i.e. the dimension of 
p is the number of control variables (i.e. the dimension of Uj^) and m 
is the number of measurements. 

In this study, the steady state gains were xased to take advantage 
of the reduction in the number of on-line computations required. If 
no advantage of the special form of a and 3 is taken, the niJiiiber of 
multiply-add operations required is 646 per update. Taking advantage of 
the form of a and 3 as given in (87b) this number can be reduced to 468. 
Further sxjbstantial reductions are possible using the properties of <|> , 

(}i , Cl and C 2 ; however these were not investigated in this study. If 
new estimates are required at a rate of ten per second, the above 
computations would require that one multiply-add operation be performed 
in a neximum period of 200 ysec, including memory access time. Sim- 
ulations of the filter under various conditions were done. Plots of 
these rjns are shown in Figures 6-12; these are discussed in Section V. 
A block diagram of the filter is given in Figure 4. 
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IV. Development of Digital Automatic Control Law 


During the final approach phase, the aircraft aligns itself 
with the runway during localizer capture and stabilizes on a strai^t 
and level course, it remains on this course until the glideslope 
capture maneuver starts, then follows the glideslope until the flare 
maneuver brings it to touchdown. Thus, as the aircraft starts to capture 
the glideslope its control surfaces are set so as to follow a constant 
altitude straight line flight path. Hence, the function of the control 
law to be developed is to take the aircraft from this condition to 
another steady condition with a constant sink rate and follow the 
fli^t path defined by the glideslope in an automatic mode using a 
digital ccmputer for the control law computations. In terms of the 
aircraft equations of motion given by (50), this corresponds to starting 
from an initial state that describes a constant speed level fli^t 
condition and bringing the state variables x to a value as close to zero 
as allowed by the wind conditions. As the state vector x represents 
the deviation of the aircraft longitudinal variables from their value 
on the glideslope, a value of zero for x means no deviation from the 
desired fli^t path. As the winds, however, will cause deviations from 
this flight path, it is necessary for the control law to take action 
against deviations caused by random wind gusts. These objectives can 
be mathematically described by a quadratic cost function that penalizes 
more for large values of x (i.e. large deviations from the glideslope) 
than values close to zero: 


38 


J(u) = ~E / [x*(t)Qx(t) + u’(t)^(t)]dt, (92) 

■^o 

where Q and R are non-negative definite matrices chosen so as to reflect 
the relative importance of deviations in Idle state variables and controls 
from their desired values and E denotes the statistical expectation 
operator. Thus, a control law is judged depending on the value of the 
cost function J(u); a small value for J(u) corresponds to small deviations 
from the glideslope, hence to a "good" control law. Thus, it is desirable 
to find a control law for which the cost function is minimum or at 
least small. 

A. Discretization of the Ctost Function 

The aircraft equations of motion (50) and the wind model (72) were 
discretized in Section IIC. The cost function (92) can also be 
discretized [15] under the assumptions that the measurement noise 
and the wind generation noise have zero mean and are gaussian, and 

that the control u(t) remains constant over each sampling period, i.e., 
equation (74) holds. If the cost function is discretized then the 
optimal control problem .of minimizing J(u) under the constraints of (80) 
and (76) can be solved using dynamic prograimiing. 

First note that (92) can be broken down to a sum of terns by in- 
tegrating over each sampling period. 

T M \+l 

J(u) = yE ^ f [x’(t)()x(t) + u' (t)Ru(t)]dt. (93) 

k=0tj^ 

Now, using (50) and (72), and x(t) and W(t) can be expressed as 
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x(tj^ + s) = (J>(s)x(t|^) + r^(s)WCtj^) + r(s)Uj^ + rij^Cs) 


(94a) 


W(tj^ + s) = (|)^(s)W(tj^) + C]^(s)j 0 < s < T, 


As 

where (J)(s) and (|) (s) are the inatrix exponentials e^ and e ^ , 
w 

respectively, and 


r(s) - [| (j)(s - T)dx]B, r^(s) = J <Ji(s - t)DC ^<J)^^(T)dT, 


w w 


n, (s) = / 


/ <J)(s - t*)DC <(» (t* - T)dx' 

•' w w 


+ x)dx, 


5]^(s) = / <P^(s - + x)dx. 

o 


Note that x(tj^), W(tj^), r(T), r^(T), <J)(T), 4>„(T) , rij^(T) and 
correspond to Xj^, W^,, r, r^, <J>, (|)^, ri^, and respectively, in 
previous notation. The above expressions for x(tj^ + s) and W(tj, 


can be substituted into a general term of the summation in (93), 
some manipulation we obtain 


\+l _ _ T _ 

Ej [x*(t)^(t) + u'(t)Ru(t)]dt = E/ [x*( 1^ + s)Qx(tj^ + s) 


u'(tj^ + s)Ru(tj^ + s)]ds 


= + 2Xj^NWj^ + 2 X|^Mu^ + 2Wj^Su^ + u^^RUj^] + , 


(94b) 

(95a) 

(95b) 

(95c) 

5 j,(T) 

+ s) 
after 

(96a) 

(96b) 
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vdiere 


Q = / ()) * (t )Q(}) (t )dt , 


(97a) 


M = J (j)'(T)Qr(T)dT, 


(97b) 


N = / 4i'(T)Qr^(T)dT, 

O 


(97c) 


s = / r’(x)Qr(T)dT, 


w 


(97d) 


R = / r’(T)Qr(T)dx + TR, 
o 


(97e) 


and is a constant depending only on the covariance of the white noise 
process C(t). It should be noted that Hj^(s) and 5j^(s) depend only on 
the values of C(t) which occur after tj^; this can be seen from (95b) and 
(95c). On the other hand and depend on the values of 5(t) before 
tj^. Since 5(t) is a zero mean gaussian white noise process, its 
values before tj^ and those after t^^ are independent; thus, 


E(Xj^nj;.(s)) = = E(^^^^(s)) = 0 


(98) 


Hence, the terms in (96a) which involve 1lie cross-correlations in (98) 
do not influence the discrete cost function (96b). On the other hand, 
note that c^ does not depend on the values of the Uj^. Thus, it does not 
affect the minimization of the discrete cost function with respect to Uj^ 
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and can be ignored. The above result is summarized in the following 
lemma. 


Lemma 1 The cost function J given by (93) and the discrete costs 
function Jj given below differ by a constant independent of u^, hence, 
J and Jj are equivalent as far as minimization with respect to Uj, is 
ooncemed. 


M 


Jl(u) = 2 


(99) 


Thus, the continuous cost function can be e><pressed in terms of the 
samples of the variables at the sampling instants, t . Note that (99) 

K1 

contains cross product terms between the aircraft variables Xj^ and the 
wind velocities and the controls u^, whereas the continuous cost 
function contains no cross product terms. It can be seen from (97) that 
the cross terms as well as the quadratic terms represent the effect of 
the dynamics in between sampling instants. As the system response 
between the sampling instants is included in the discrete cost function, 
it becomes possible to use sampling periods larger than generally used, 
thus allowing more time for computation during updates of the estimates 
and the control. 

B. Solution of the Optimal Control Problem 

The optimal control problem of minimizing a quadratic cost function 
with the constraint of a linear dynamical system with gaussian statistics 
has been extensively treated in the literature [22 and the references 
therein]. The problem considered here differs from the usual one in two 
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ways; the cost function contains cross products terms as mentioned in 
the previous section and the state of the system is affected by a 
random disturbance which can be described by linear dynamics; i.e. the 
disturbance is not white noise. The case of non-random disturbances 
was treated in [23], [24], [25]. The problem considered here is also 
treated in [26]. 

Consider the optimal sampled-data control problem posed by the 
aircraft dynamics (50), the wind model (72), the control constraint (74), 
the discrete measurement equations (85) and the costs function (93). In 
earlier sections this sampled-data problem has been reduced to a 
discrete form with the discrete aircraft equations described by ( 80 ) , 
the wind model by (76), the measurements by (85) and the cost function 
by (99). Thus, the optimal control problem is reduced to minimizing 
the cost function (99) with respect to the control sequence where 

Uj^ is restricted to depend only on the past measurements {y^, 0<i<k-l} 
as these measurements only are available for controlling the aircraft.* 

To derive the control Uj^ which minimizes the cost function (99) the 
dynamic programming method [27] will be used. Thus, consider a general 
term in (99) . 

= i. * 'Hh? 

( 100 ) 


* This restriction is nathemtically interpreted as: ^ is measurable 

with respect to the a-algebra generated by the random variables 
{y. , 0 < i < k - 1}. ' 
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Note that as E(x*Qx^) and E(x*NW ) are constants they can be excluded 
from the minimi 2 ation, so that the general term in (99) can be written 
as shown in (100). Now, substituting (80) and (76) into (100) and 
regrouping terms 


= I * Hi ♦ '“k”k 

* KW * 


( 101 ) 


vdiere 


D, n P,, r + P„, 4 > , R, n R + r'P-,, T, 
k Ik w 2k^w’ k . Ik ’ 


(102a) 


®ik = ’"Pik* " Sk “ '•'“k * ' 


(102b) 




(102c) 


First note that the cross-correlations of x^^, ^ with Dj^ or 


are zero as shown in (98), and 
k ’ 


E(u^n') = E(E(u^,’|/^» = E(u^E(,;|>-^)) = 0< 


(103) 


and similarly for E (u^ ; so that these terms have been dropped from 

(101). Also note that d^, depends only on the statistics of the 

* y, is the a -algebra generated by the past measurements {y., 0<i<k-l), 

m. ^ 
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disturbances (i.e. winds), but does not vary witii the control sequence. 
Now, note that since Uj^ depends only on past measurements, using a 
well-known lemma on conditional expectations C8], 




^ /N 

where Xj^ and Wj, are the conditional expectations ECx^li/j^), E(W^|Vj^) 
respectively. Note that Xj^ and are the components of ^ in (89) so 

A. A. 

that Xj^ and are obtained from the filter output. Substituting (104) 
into (101) , 


= f * A* '“A * * SA> 

(105) 

Minimization of with respect to u^, results in 




(106) 


Note that since x^^ and depend only on the past measurement, so 
does and the restriction on Uj^ is thus satisfied by (106). Now, 
substituting (106) into (105) and adding the next general term in (99), 
we obtain 
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^k-1^^-1^ " 2 ^^-l^uk-1 ^^k-l^\-l 

(107) 

^ K-iK-i^ ^ Vr 


where 


'’ik-l “ ♦■'’uc* - ®ij<®k‘®uc + ‘5- f'M = 


(108a) 


"’2k-l = ^2M = N 


(108b) 




(109a) 





(109b) 


Since, Wj^ and are gaussian random variables, it is known that 
[21] E(X|^x^), E(Xj^W^) and E(Wj^W^) do not depend on the value of u^; hence, 
^ is constant as far as minimization with respect to u^. Since, 

(107) is of the same form as (100) except for a constant, the above 
steps can be applied to and continue the iteration until is 

obtained. Hence, the following theorem summarized the result obtained. 


Theorem 1 The optimal control problem posed by (50), (72), (74), (85) 
and (93) has a unique solution (a.e.) which is given by (106), (108) and 
( 102 ), 
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( 110 ) 


"k = = -«lk^ - ^2k\* 


The control is se^ to consist of two parts. One uses the aircraft 
state estinates for feedback; note that tiie gain associated with 1diis term 
is the same gain as the case when no disturbance is present. Thus this 
term accomplishes the stability about the glideslope to be followed; 
i.e., if the aircraft is not on the glideslope, its effect will be to 
bring the aircraft back on the glideslope. The second term uses the 
disturbance estimates; it acts as a preventive measure so that if a 
disturbance that will throw the aircraft off the glideslope is estimated 
t>y corrective control action will be taken before the aircraft is off 
the glideslope, and thus prevent fli^t path effors as much as possible 
before "Uiey occur. 

C. Constant Feedback Gains 


From the expression (11) for the control, it is seen that the feedback 
gains and vary with time, i.e. with k. Implementation of this 
control law would require solving the equations (108), (102) backwards, 
starting at K n M, to determine the gains at each sampling instant, tben 
storing these gains on an on-board computer to use them in the computation 
of the control. This would place a great burden on the computer's 
memory capacity and speed. The use of these time varying gains, however, 
is not necessary to obtain a good performance and stability. It is 


possible to use constant values for the feedback gains and thus 

avoiding implementation conplexity. In fact, with the use of constant 


gains to determine the control in (11) , the control ccaiputations require 
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very little complexity in the on board computer. In this case, it is 
necessary to store idle elements of and and compute (11) within 
a time period T, i.e. the sanpling period. For the system considered 
here, the number of elements in and H 2 is 45, and the number of 
multiply-add operations required for each update of the control is 45. 
Thus, the control computations can be easily implemented on an on board 
computer. 

To compute the constant gains to be used, it is necessary to solve 
equations (108) and (102) backward. It is known that converges to 
a steady value under very loose conditions [21], which are satisfied by 
the aircraft equations ( 50 ) . The convergence of P 2 J, when ())^ has unstable 
poles is not always guaranteed. The following theorem given here without 
proof specifies necessary and sufficient conditions for the convergence 
of P 2 j^. Let p((j)) denote the largest of the absolute values of the 
eigenvalues of cj). 

Theorem 2 Let P^,_ converge. Then P„,_ given by (108) converges if, 

-LK 2k 

and only if, 

p((j) - rH^)p(ij)^) < 1 

Thus, the convergence of P 2 ^^ and H 2 k depends on the degree of 
instability in the disturbance model. For the wind model considered here 
p(<j)^) has a value of 1, v^ile p(<t> - FH^) is less than 1. Hence, the 
above condition is satisfied, and P 2 Jc’ ^2k converge. The steady 

value of P 2 j^ and the gain H 2 j^ can be conputed from (108) and (102). 
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The steady state value of was conputed using a non-iterative 
approach [28] ; i.e. , the steady state value, say P^, vias con 5 >uted without 
using (108a) vAiidh requires a large number of iterations to converge. 

In any case, this computation would be done off-line and place no real- 
time requirements. The steady state value of P^j^ was also computed using 
a non-iterative method, developed by the author, which will not be given 
here. This method also reduced off-line computation time. Then, the 
gains and H 2 were determined using (11) and (102b). These gains were 
stored and used in the simulation of the control system. To avoid 
extremely high control coimiands which mi^t be caused by component or 
sensor failures and would place hi^ stress on the control siirfaces 
limiters were placed as shown in Figure 5. Note that the stabilizer and 
throttle commands (u 2 j^» u^j,) are rate coimiands. If the stabilizer and 
throttle positions are needed, a hold circuit followed by an integrator 
can be used. 
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V. Results 


Using the models developed in the earlier sections a versatile 
digital computer simulation was developed. The filter and the control 
law developed for the capture of a steep 6° glideslope were simulated on a 
digital con 5 >uter. The aircraft equations of motion, linearized about a 
6° glideslope developed in Section II. A and discretized for simulation 
in Section II. C , were coded on a digital computer to simulate the motion 
of the aircraft. The wind model developed in Section II. B was used to 
generate random gusts as well as steady winds to simulate the wind 
conditions for the simulation of a given flight. The MLS receiver outputs 
and on board sensor outputs were simulated by corrupting the position, 
velocity and attitude of the aircraft by noises characteristic of present 
day sensor errors, the models for these sensor errors are given in 
Section III. A. Thus, the simulation developed included the following: 

1. aircraft motion in the longitudinal axis 

2. wind conditions (gusts and steady winds) 

3. sensor errors 

4 . filter 

5 . control law 

Major parameters such as turbulence levels, steady winds were left 
as input variables in the simulation developed to allow for versatility 
of use. Thus, to simulate the aircraft motion under different wind 
conditions it is only necessary to specify the values of the wind gust 
and steady wind velocity parameters. The mjor parameters that can 
be specified include: 
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1 . 


the standard deviations of wind gusts 

2. scale of turbulence 

3. steady wind velocities 

4. standard deviations of the noises of each sensor and MLS 

5. initial conditions at the start of glideslope capture 

6. sampling rate, T 

Note that by varying the sampling rate, T, it is possible to simulate 
varying rates for the reception of MLS data; by varying the noises on the 
various sensors, it is possible to degrade the quality of the MLS signal 
by adding noise to it and similarly for the other sensors. 

A. constant gain Kalman filter was developed to filter out the noise 
present in the various sensor outputs and to obtain estimates of various 
unmeasured parameters, such as wind velocities, and angle of attack. 

The filter combines position data from the MLS with air data from 
on-board sensors, body-mounted accelerometer data and aircraft attitude 
for this purpose. The estimates output by the filter are used to canpute 
the control commands. The simulation runs made under varying wind 
conditions indicate that- the estiiiates supplied by this filter follow 
the actual aircraft parameters with little error. 

A digital automatic control law for the longitudinal axis to 
capture and follow a 6° glideslope was developed. The design procedure 
used was the general quadratic regulator theory with the modifications 
described in Chapter IV. The control law uses constant gains in the 
feedback loop, and also uses wind velocity estinetes provided by the 
filter to retrim the aircraft and reduce the deviations caused by the 
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wind gusts. The 4D aspect of the control law is due to the fact that 
the position error fed back is the difference between the actual aircraft 
position and its desired position at each instant of time. Thus, the 
■aircraft control law tries to bring the aircraft to a specified position 
at a specified time. The control law was simulated under various wind 
conditions and from different initial conditions as shown in Figures 6-12. 
The 4D errors in position, errors in sink rate and inertial speed are 
shown in the various plots. The performance of the aircraft in captxrring 
a steep 6° glideslope from level fli^t at 120 knots appears to be 
acceptable. 

Various simulation runs for capture and glideslope following are 
shown in Figures 6-12. The initial conditions corresponding to time 
zero, 30,000 ft. from the runway and an altitude of about 3,153 ft. have 
been specified as the beginning of the capture mode; i.e. at time zero 
the aircraft is trimmed for level flight at 120 knots (202.536 ft/sec); 
the capture mode is switched on when the aircraft's estimated altitude 
below the glideslope is 45 ft. Thus, the initial 4D vertical error 
shown in the simulation plots is not interpreted as an error but as a 
starting point from which the glideslope is to be acquired. Similarly, 
the initial conditions for pitch, sink rate, thrust, stabilizer, etc. are 
at the values required for level fli^t. The units used in the plots 
are ft. for distances, ft/sec. for velocities, degrees for angles, 
degrees per second for angular rates and seconds for time. Thrust is 
expressed in pounds, throttle in degrees (on the stick), elevator, and 
stabilizer in degrees. During the initial capture period, the 
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aircaraft acquires the 6° glideslope by retrimming the stabilizer position, 
vAiile the thrust is reduced. The elevator initiates the pitch-down 
notion, necessary to acquire a sink rate (h) of about 21 ft/sec. The 
settling time for pitch, pitch rate, sink rate is about 10 seconds, 
while the position errors and speed take longer to reach their steady 
values. No extreme values are required for the controls during capture, 
and overshoots are not excessive .althou^, in comparison to 1±ie usual 
3® capture, the 6° capture is a relatively large change in fli^t 
condition. 

To compare the effects of various levels of wind gusts, steady winds 
and sensor errors, tiie hypothetical case of no sensor noise and no wind 
was simulated (Figure 6). The aircraft notion for tbis case is smootb, 
and after capture tbe deviations from the glideslope settle to zero. The 
introduction of sensor noise on the measurements (Fig. 7) causes sli^t 
but noticeable deviations from the glideslope in all the variables; 
this is caused throu^ the errors in the estimates of the aircraft 
variables due to imperfect measurments. The effects of wind gusts on the 
aircraft response can be seen in Figures 8-10. The turbulence levels 
are specified throu^ the standard deviation of longitudinal gust 
velocities, a , and the vertical gust velocities, o . As can be seen 
from the plots, the wind gusts affect the aircraft motion considerably. 

The control action increases in order to reduce the fli^t path deviations 
from the glideslope by using the proper feedback. The maximum deviations 
while following the glideslope are within acceptable limits in the 
various wind conditions. It shoiold be noted that these deviations could 
be further reduced by allowing a hi^er level of control activity. 
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This can be achieved simply by reducing the value of the elements of R 
in the cost function (92). To see the effect of a hi^er sanpling rate 
on the performance of the aircraft motion, a rate of 10 per second 
instead of 5 per second was used for incoming data and control commands. 
Ihis case is shown in Fig. 12; a sli^tly titter performance in the 
aircraft deviations can be seen; however, the basic response character- 
istics are the same. It should be noted, however, that even though 
rates of 5 per second appear to be satisfactory for glideslope capture, 
the flare maneuver may require a hi^er rate. In general, the control 
law appears to be performing satisfactorily during the glideslope 
capture and glideslope following phases of final approach. 

The various plots with sensor noise and wind gusts show that the 
constant gain filter tracks the various parameters with considerable 
accuracy. The case for no wind and no measurement error (Fig. 6) is 
given as a reference for comparison of the effects of inducing sensor 
noise and winds of various mgnitudes. The case of no measurement noise 
(Fig. 6) shows no noticeable error in the estimates, as would be 
expected. When measurement noise is introduced (Fig. 7) errors in the 
estimates become noticeable. When wind gusts are introduced (Fig. 8), 
the errors in the estimates increase sli^tly; however, this increase in 
the errors is small in comparison to the errors due to sensor noise 
(Fig. 7). Also note that the wind velocities are estimated with 
accuracy, the measurements containing information about wind velocities 
are airspeed and the accelerations; these appear to be sufficient to 
track the wind velocities. Steady winds are also estimated with accuracy 
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as can be seen in Fig. 9. Figure 11 shows a case where the initial 
values of the estimates are different than the actual values. In 
general, it can be said that the constant gain Kalman filter has a 
satisfactory performance in filtering sensor noises out and providing 
estimates of unmeasured variables. 
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VI. Conclusions 

In the previous sections a digital automatic control law for steep 
glideslope capture and glideslope tracking was developed. The control 
law consists of a filter and a controller. The filter accepts MLS data, 
air data, attitude data and body-mounted accelerometer data, filters 
out the sensor noise in the measurements and outputs estimtes of the 
measijred variables, aircraft velocities and wind velocities. The 
controller uses these estimates to compute surface setting commands. 

A digital computer simulation of the aircraft motion, sensor noises, 
wind conditions, the filter and the controller has been developed to 
test the concepts used in the development of the overall control law 
for automatic steep glideslope capture and tracking. On the basis 
of the simulation results, the use of steep glideslopes during the 
glideslope capture and glideslope tracking phases of the final approach 
appears to be feasible. 

For the filter concept used in this study, the use of costly 
Inertial platform data can be replaced by less costly body-mounted 
accelerometer data (which have less accuracy) to obtain sufficiently . 
accurate estimates of the variables needed for control purposes. It may 
be possible to eliminate the use of accelerometer data using this type 
filter as long as accurate MLS data is available; however, as inexpensive 
accelerometers are available, this case was not considered here. 

Further reseaipch is required to determine if the accuracy of velocity 
and acceleration estimates obtained using only MLS data would be 
sufficient for successful autonatic landing under turbulence. The 
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simulation results also indicate that the concepts used in the develop- 
ment of the control law can provide successful control action for steep 
glideslope capture and tracking under various levels of turbulence. 

Another critical phase of the final approach to landing is the 
flare maneuver. During a steep approach, due to larger changes in pitch 
and sink rate required, the flare naneuver is more critical than for 
shallower approaches. Thus, further research to investigate the problems 
which nay be encountered in flare during a steep approach would be 
required for a complete evaluation of steep approaches to landing. 
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Figure 5 Block Diagram of Control System 

(c is a discrete specifying capture; 
c means the control system is not 
in the glideslope capture mode) 
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Table 1 Standard Deviation Values for the 
Simulation of Sensor Noises 


Variable 

Standard Deviation 

Type 

Yi 

.15° 

additive 


.10° /sec 

additive 

Ys 

1 ft. 

additive 


.031° 

additive 

Ys 

25 ft. 

additive 

Ye 

5 % 

multiplicative 

Yy 

.005 g 

additive 

Ys 

2 % 

multiplicative 

Yg 

.005 g 

additive 
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Appendix A 

The aircraft equations of motion vdiich are used in the simulation 
and for the development of the filter and control law were developed in 
Section II. A. The final form of the continuous equations as given in 
equation (50) is repeated here for convenience 

X = Ax + Bu + rw 

The form of the matrices A, B and D are also given in Section II .A. 
In this appendix, we shall give the expressions for the non-zero elements 
of these matrices in terms of the aircraft stability derivatives; the 
stability derivatives are assimied to be in the stability axis. 

mg COS0 mg sin0 

aj = i CX2 - ~ > 

03 = , 04 = <03 + C^) , 

“5 = = I P^O 

where S is the wing area, c the mean aerod 3 mamic chord and p is the air 
density. Using these variables, the matrix elements are given below: 
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aj 4 - 1 


” ao' ’ ^22 = ~ Sboi * ^^Txo^ ’ 


^23 


CI 3 as 

"lo 




«3 


-, _ — . - D6s 

’ ^27 - — , a 28 - -— , 


^31 ~ “ 2*^4 » ^32 “ '*' ~ '*' ’ 


La Do 


I 34 = a 4 (o 3 - Cj^) , a 37 = « 4 C ^2 ’ %9 ~ “‘**"L 6 s ’ 


a 20 t 4 -] 

3 . 4 -I — — C • 3 ^U 9 “ (C 2C ^32^ * ^ a 

as ma as mu mo ma 


343 - (C + C rp ^3 3^ * » ^44 “ /-, ^ ^34^ * ^ 5 

as im iriTa ma ’ as mq ma 


a 47 = 


C rp ■*■ ^376 • 

nfT ^ ma 
«5 


^49 = 


C . + &4 qC • 

m6s mg 

«5 


a 7 7 — — . 5 , 3.J 0 - .298 j 


^ _ *"D 6 e T- _ o K - ^m 6 e ^ 

^21 ^ » 1^31 - -“4CL5e ’ ^41 


^33 = « 4 (C^ - > d 43 = Jia. 


C + (dqo - DC 


33 


ma 


Thus, given the stabDity derivatives of the aircraft, the A, 
D natrices can be computed. 


B and 
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